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Abstract 

We consider the effect of radial fluid injection and suction on Taylor-Couette flow. 
Injection at the outer cylinder and suction at the inner cylinder generally results in 
a linearly unstable steady spiralling flow, even for cylindrical shears that are linearly 
stable in the absence of a radial flux. We study nonlinear aspects of the unstable 
motions with the energy stability method. Our results, though specialized, may 
have implications for drag reduction by suction, accretion in astrophysical disks, 
and perhaps even in the flow in the earth's polar vortex. 
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I. THE PROBLEM 



There are many contexts in which rotating flows arise and their instabilities are naturally 
of interest. For an axisymmetric, incompressible, inviscid flow with no variation along the 
axis of symmetry, Rayleigh (1916) provided a necessary and sufficient condition for stability 
against axisymmetric perturbations (Chandrasekhar, 1961). Rayleigh also showed that a 
necessary criterion for the instability of such flows to non- axisymmetric perturbations is an 
analogue to his classic inflection point criterion for plane flows. He expressed this criterion 
in words but the correct mathematical expression is easily worked out following Rayleigh's 
lead (Spiegel and Zahn, 1970). For the study of the effect of viscosity on these results, 
the Taylor-Couette configuration of flow between coaxial cylinders provides a good proving 
ground (Chandrasekhar, 1961), although shear and rotation come together in other contexts 
as well (Yecko and Rossi, 2004 and references therein). 

When criteria for the simplest two-dimensional circular shear flows indicate stability, the 
question naturally arises whether some simple extraneous effects may be destabilizing. In 
the case of the linearly stable plane shear flow, magnetic fields (Stern, 1963 and Chen and 
Morrison, 1991), Rossby waves (Lovelace et al, 1999) and suction (Hocking, 1974; Doering 
et al, 2000) have been found to be destabilizing mechanisms, given suitable conditions. In 
the case of a stable Taylor-Couette flow, it had been seen even earlier that a magnetic field 
may be destabilizing (e.g. Chandrasekhar, 1961; Balbus and Hawley, 1998). In an attempt 
to add to this compendium of instabilities, we here investigate the instability that arises 
when a formerly circular flow is turned into a spiralling flow by the imposition of continous 
suction on the inner cylindrical boundary with a compensatory injection of fluid on the outer 
boundary. For brevity, we shall hereinafter refer to the induced radial inflow as accretion 
and the mechanism driving it as suction. 

A practical example of the conversion of differential circular flow to spiraling motion, or 
swirl, arises in the study of drag reduction on a cylindrical airfoil (Batchelor, 1967). In that 
case, when suction is imposed on the central cylindrical airfoil, the drag is reduced because 
boundary layer separation is inhibited. Spiralling flows are also central to the structure and 
dynamics of accretion disks that form in a great variety of astrophysical situations such as 
binary stars (which are often X-ray sources), flows around massive black holes at the centers 
of galaxies, and in the disks around newly born stars where planetary systems are thought 
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to form. 

Of course, accretion disks are more complicated than the simple model that we consider 
here. In most astrophysical cases there is ionization so hydromagnetic effects become im- 
portant and can produce the magneto-rotational instability (see Balbus and Hawley, 1998). 
However, protoplanetary nebulae such as the primitive solar nebula are believed to have 
been cool so that other mechanisms may need to be considered to get the full picture (for 
at least some phases) of the turbulent dynamics of planet formation. The feature of the 
present study that may have relevance to the fluid dynamics of astrophysical disks is the 
radial inflow, or accretion that, in some astrophysical examples, is fed by inflow from an 
external mass source. A similar combination of motions occurs also in the terrestrial polar 
vortex that may play a role in the transport of ozone in the upper atmosphere (Mclntyre, 
1995). 

With these motivational examples in mind, we turn to an idealized pilot problem to learn 
what may happen when rotation and accretion combine to produce a spiralling motion as in 
Taylor-Couette flow with suction. Although we cannot draw immediate conclusions about 
the applications of our results to the motivational examples just described, especially to 
accretion disks, we regard it as suggestive that the combination of rotation and accretion 
(i.e., radial injection and suction) in a flow leads to instability. As we shall see in what 
follows, taking into account both accretion and rotation leads to very different stability 
properties of the flow than when only one of these effects alone is included. For example, 
in the classical problem of the linear and nonlinear stability of Taylor-Couette flow between 
concentric cylinders, the system is linearly stable when the outer cylinder is rotating and the 
inner one is stationary. Likewise, simple radial accretion without rotation is linearly stable. 
However, we find that, at sufficiently large rotational Reynolds numbers in a stable TC flow, 
small accretion rates lead to linear instability. Moreover we find that the rotational motion, 
rather than simple plane shear, plays a major role in the stability characteristics. 

An initially plane-parallel shear layer perturbed by transverse suction may also exhibit a 
small-suction instability (Hocking, 1974), while larger suction rates absolutely stabilize the 
corresponding laminar flow for arbitrarily large Reynolds numbers (Doering et al, 2000). 
This strong stabilizing effect of suction for even very high Reynolds numbers in the plane 
case does not carry over to the rotational case; in the cylindrical geometry transient growth 
of perturbations is generally encouraged by strong suction. 
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The following discussion of these issues is organized as follows. In section [TT] we describe 
Taylor-Couette flow with suction. Then we treat the linear stability of this flow in section 
II III where we find linear instabilities in the simple flow of the model problem. In section 
IIVI the asymptotic behavior of the onset of linear instability is computed for different limits 
in the narrow-gap approximation. We next go on to consider some nonlinear aspects of 
the instability. First, in section |V] we describe the energy stability method and then derive 
some rigorous bounds on the energy stability limit of the flow in section IVII That limit is 
computed numerically in section IVII I The concluding section IVIIII summarizes the results 
and mentions some possible implications. In an appendix we note that physically relevant 
bounds on the mean energy dissipation rate for non-steady flows (including statistically 
stationary turbulence) remain unknown in the Taylor-Couette geometry with suction. 

II. FORMULATION OF THE PROBLEM 

A. Geometry of the problem and boundary conditions 

Consider the incompressible flow of a Newtonian fluid between two coaxial porous cylin- 
ders. For definiteness, we let the inner cylinder be at rest while the outer cylinder is rotating 
with constant angular velocity, Q. The radii of the inner and outer cylinders are, respec- 
tively, Ri and R2, and we use cylindrical coordinates (f, 9, z). We consider periodic boundary 
conditions in z with period L z , and impose no-slip boundary conditions so that the vertical 
(z) component of the velocity field vanishes at both cylinders and the azimuthal component 
matches the velocity of each of the cylinders. 

Now we add a radial inward flow to this classical Taylor-Couette configuration. Fluid is 
injected at the outer cylinder with the constant volume flux 2inp (volume of fluid injected 
per unit time and per unit height of cylinder) and uniformly sucked out at the inner cylinder 
at the same rate. The control parameter (p > is central in what follows. The resulting 
velocity field is u = ue r + veg + we z and the boundary conditions are 



u(R 2 , 0, z) 



u(i?i, 9, z) 



= -TT e r + R2^Gg. 
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(1) 



(2) 
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The incompressibility constraint V.u = implies two further boundary conditions: 

u f (Ri) = and u f (R 2 ) = (3) 

where the subscript f denotes partial differentiation with respect to r. A schematic of the 
set-up is shown in Figure [TJ 




FIG. 1 The outer cylinder rotates at angular velocity £1. Fluid is injected at the outer boundary 
with an entry angle O = arctan^/di?^)], and removed uniformly on the surface of the inner 
cylinder. 

B. Dimensionless numbers 

This problem involves five parameters: the two radii, R\ and R2, the angular velocity 
of the outer cylinder, Q, the imposed accretion flux, 2inp, and the viscosity of the fluid, v. 
Nondimensionalizing time and space leaves three independent dimensionless numbers rj, G 
and Re describing the geometry and physical boundary conditions: 

• The geometrical factor is 77 = When 77 — > 1 we approach the narrow gap 

limit where (R2 — Ri) << Ri, and expect to find results similar to those from 
the slab geometry — namely plane Couette flow with suction. On the other hand, 
when 7] goes to zero, the outer radius R2 goes to infinity relative to the inner ra- 
dius and we expect to see significant effects of the curvature on the stability of the flow. 
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The accretion flux is measured by the injection angle, 6, at the outer cylinder, where 

tane = 7§r < 4 > 

If tan© = 0, there is no suction and we recover a classical Taylor- Couette problem 
with a linearly stable stationary flow when the inner cylinder is stationary. When 
tan 6 — > oo we approach the limit where suction dominates rotation. 



The azimuthal Reynolds number is 

Re = R^R.-R,) = R 2 2 n(l- V ) 

V V 

We choose this definition of the Reynolds number to match the definition of the 
Reynolds number for the plane Couette flow when 77 — >■ 1. 



• Introduction of an additional parameter, the radial Reynolds number a = ip/v, is use- 
ful for making some of the equations more compact. This parameter is not independent 
but is linked to Re, i] and tan(O) by the relation 

i?e tan 6 . . 

We define the dimensionless velocity u = ue r + ve$ + we z = u/(R 2 £l) and the dimension- 
less coordinates r = f/R 2 and z = z/R 2 and the time t — Qt . The velocity field satisfies 
the incompressible Navier-Stokes equations, written in dimensionless form as 

u t + u-Vu = -Vp + ^Au (7) 

Re 



with boundary conditions 



, . . tan0 

u(r],9,z,t) = e r (8) 

V 

u(l, 8, z,t) = — tan6e r + eg (9) 

u r (r],e,z,t) = (10) 

u r {l,9,z,t) = (11) 



6 



C. The basic laminar solution 



The basic steady laminar solution Yi am = (U{r), V(r), 0) of equation ([7]) has radial 
velocity 



TT , . tan 6 , . 

u (r) = — — (!2) 



and azimuthal velocity 



V(r) = ^r l ~ a + . (13) 

The azimuthal velocity profile is represented in Figure [2] for several different combinations 
of values of a and rj. 

When a — > we recover the classical velocity field of the Taylor-Couette flow with a 
stationary inner cylinder. That flow is independent of the kinematic viscosity, u, and the 
azimuthal velocity then increases monotonically when r increases from rj to 1. 

For nonzero values of the entry angle tan and for sufficiently large values of the Reynolds 
number, that is, for .Re tan O >> 1, the azimuthal speed increases outward (from the inner 
boundary) in a boundary layer of thickness 5. On the other hand, as the fluid comes in from 
the outer boundary, its motion is mostly inviscid and thus conserves angular momentum 
rV{r): the azimuthal velocity of the fluid element increases during the inward journey from 
the outer boundary according to 

V(r) « - (14) 
r 

For a > 1, the two flows join to form an azimuthal velocity profile with a maximum at 
r = 7] + 5: 

V r (r] + 5) = 5 = r][(a- 1)^ - 1]. (15) 

This particular feature of the azimuthal velocity profile, that it can achieve a local maximum 
within the gap for a > 1, is fundamentally linked to the rotational geometry. This effect 
disappears in the plane limit rj — > 1 since the maximum value of the azimuthal velocity 
inside the gap V max = -. For large values of a, the boundary layer is very thin; 5 ~ f] 1 -^^-. 

These considerations point up a significant but not widely appreciated feature of accre- 
tion flows. The asymptotic velocity profile as Re — > oo is dramatically modified from its 



structure at tan0 = when tan 6 7^ 0. When tan0 = 0, the azimuthal velocity profile is 
the Reynolds-number-independent Taylor- Couette profile. When tanO 7^ the azimuthal 
velocity profile varies as - everywhere inside the gap and outside the boundary layer, which 
becomes infinitely thin in the limit Re — > 00. This shows that even if the accretion rate 
is very low (but not zero), one cannot neglect it at high Re and study the problem with 
rotation only. The limits Re — > 00 and tan — > do not commute. 




FIG. 2 Azimuthal velocity profiles for different values of the radial Reynolds number (top : n = 0.9, 
bottom : n = 0.5). 
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III. LINEAR STABILITY ANALYSIS AND NON-AXISYMMETRIC DISTUR- 
BANCES 



To study the destabilizing effect of the inward suction on the stability properties of Taylor- 
Couette flow, we perform a linear stability analysis. It is already known from work by Min 
and Lueptow (1994) that, for axisymmetric disturbances to Taylor-Couette flow, suction is 
stabilizing and so increases the critical Taylor number at which Taylor vortices first appear. 
Therefore we here consider the linear stability theory for non-axisymmetric perturbations. 
We continue to keep the inner cylinder stationary, so that the flow without suction is linearly 
stable according to Rayleigh's criteria for inviscid perturbations. 

In most instability mechanisms, only one symmetry of the initial problem is broken at the 
onset of the primary instability while others are broken through secondary instabilities. Since 
we anticipate instability to nonaxisymmetric perturbations, we shall not, in the first instance, 
break the invariance to translations in the ^-direction in seeking the primary instabilities of 
the spiralling flow. We observe that, consistently with this supposition, a Squire's Theorem 
exists in the plane-shear-with-suction situation, that is, in the limit r/ — 1 (Doering et ai, 



A. Linearization of the equations 

We make the decomposition u = V iam + v for the full velocity field, write the r, 9, and 
z components of a mode of the perturbation as 



and seek the (possibly complex) eigenvalue A. The linearized versions of the Navier-Stokes 
equations for the r and 6 components of v for this mode are 



2000). 




(16) 



Am 



1—7] 

Re 

1-77 



tanB + ^2 



U r + AiU + Z]V + p r 



(17) 



u 



rr 



r 



Xv = 



tan0+ l_2Z 



v r + A 2 v + Z 2 u + 



imp 



(18) 



Re 



■v. 



rr 



r 



r 
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where p is the pressure perturbation and 



V tan 6 1 — 77 m 2 + 1 
A\ = im 1 h 



Re r 2 



V tan 6 1 — n m 2 + 1 
A 2 = im h 



Re r 2 



(l-7])2im 2V 
Zi = — , and 

Re r z r 



V r + - 



V (1 — rf) 2im 



Re r 



2 ' 



(19) 
(20) 
(21) 
(22) 



The continuity equation becomes 



imv + (ru) r = 0, 



(23) 



providing an expression for v that may be inserted into ( 1181) to give p in terms of u. Then, 
upon differentiating p with respect to r and inserting the result into (TT7|) . we arrive at a 
fourth order ODE for u(r): 



X [—r 2 u rr — 3ru r + {m 2 — l)u] 

(1-77) 



+ 



+ 



Re 



+ tan G r 



Re 
Re 



U r 



(6-m 2 ) + 3tan9-r 2 A 2 



u. 



tan9 + ^ 2 2/ , x / 
^—m + imrZi + imrZ2 — r (A 2 ) r — 3rA 2 



■it,. 



+ [m 2 Ai + imZ\ + 7777Z2 + imr{Z<2) r — r[A%) r — Ai\ u. 



(24) 



The boundary conditions are that u and 77 r vanish at r = 77 and r = 1. 
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FIG. 3 Top: Linear stability limit of several modes for rj = 0.9 (left) and angular velocity of 
propagation of the unstable mode at onset (right). Thick dashed line : m = 1; dashed line : 
m = 3; dash-dotted line : m = 10; solid line : m = 30; dotted line : m = 100. The energy 
stability limit has been added to make the comparison easier (thick solid line). Bottom: Contour 
plot of the streamfunction of the unstable mode at the onset of instability for m = 10, r] = 0.9, 
and tan© = 0.001 (left) and tan0 = 0.01 (right). The smaller the injection angle, the more the 
velocity gradients are pressed up against the injection boundary. 
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B. Numerical computation of the linear stability limit 

Equation fl2^|) presents an eigenvalue problem which can be solved numerically using a 
standard finite difference method. The linear stability limit of different modes is represented 
in Figures E] and S] for 77 = 0.9 and r\ = 0.5. There is linear instability even at small injection 
angles. The modes that become unstable for the lowest critical Reynolds number have 
m ~ j^-, corresponding to closed and almost circular cells (in the (r, 6) plane) between the 
two cylinders. As expected for a system that lacks reflection symmetry in the ^-direction, 
the unstable mode is a travelling wave that propagates in the azimuthal direction. The 
angular velocity of propagation Im(A)/m at the onset of instability is plotted in Figures [3] 
andU 

At a certain value of the injection angle, the linear stability boundary has a vertical 
asymptote in the tanB — Re plane (see Figure [3]). Above this critical value of the injection 
angle, the flow is linearly stable at any Reynolds number. Comparing the linear stability 
limits obtained for the two values of 77, we notice that the flow is more unstable linearly for 
the geometry with the highest curvature. From rj = 0.9 to 77 = 0.5, the minimum value of 
the Reynolds number for marginal stability — the critical Reynolds number — goes from 
approximately 10 4 to 2 x 10 3 and the maximum value of the injection angle that allows a 
linear instability increases by almost an order of magnitude. A few plots of the unstable 
modes at the onset of instability are presented in Figures [3] and HI Close to the critical 
injection angle the cells occupy the full width of the gap, whereas for smaller values of 
the injection angle the modes are pressed up against the injection boundary. In this latter 
situation, the angular velocity of propagation tends to 1, corresponding to the unstable mode 
being advected by the flow at the dimensionful angular velocity Q of the outer boundary. 

As already remarked, these linear instabilities exist when both rotation and accretion are 
present but they are lost if only one of these two ingredients is present in the configuration 
considered here. 

IV. NARROW-GAP LIMIT 

Since the numerical problem is readily solved in the linear case, we have the luxury of 
using some rough and ready approximation methods to extract analytic results about linear 
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FIG. 4 Top: Linear stability limit of several modes for rj = 0.5 (left) and angular velocity 
of propagation of the unstable mode at onset (right). In this high curvature geometry, linear 
instability is found for lower values of the Reynolds number and up to higher values of the injection 
angle. Solid line: m = 1; dashed line: m = 2; dashed-dotted line: m = 3; dotted line: m = 4; thick 
solid line: energy stability limit. Bottom left: Contour plot of the streamfunction of the unstable 
mode at the onset of instability for m = 2, r/ = 0.5, and tanO = 0.005. The streamlines are 
concentrated toward the injection boundary. Bottom right: Contour plot of the streamfunction 
of the unstable mode at the onset of instability for m = 2, rj = 0.5 and tanG = 0.1. The cells 
occupy the whole gap. 
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stability on which one may build intuition about the solutions. The numerical results provide 
guidance in the choice of effective approximations and the approximate solutions obtained 
serve to confirm those results. We next carry out such approximations in the narrow-gap 
limit. We consider situations where the critical Reynolds numbers are high enough that the 
laminar velocity profiles can be approximated throughout the gap by the their asymptotic 
forms for high Reynolds number, 

— tan G/r 

1/r . (25) 


This solution does not have the thin boundary layer of the large Reynolds number case. 
Hence we do not satisfy the no-slip boundary condition on the azimuthal velocity at the 
inner boundary in this approximation. 

The vorticity, u, of a weak two-dimensional disturbance to this basic flow obeys the 
advection-diffusion equation, 

d t u + (V a .V)u; = ^Aw. (26) 
tie 



The velocity perturbations ffl6|) imply that the vorticity perturbation is of the form u = 
Cj{r)e~ xt e imd . Introduction of this expression leads to 

1 — 77 / „ 1 — ri\ u r ( im m 2 1 — ri \ , . 

_^__i^^ tm0+ _i^ + ^ + __iy =o . (27) 

We introduce the variable x = 1 ~ r _^ 2 G [0, 1] and define parameters S, e and ii according to 

\2 



S = m(l-r 1 yRe (28) 

9 tan0 , . 
m(l — r\Y 

A = im + fim(l — rj). (30) 



Equation ( 1271) is then equivalent to: 



M^)t^ (31) 

/ m 2 (l - ri) 2 M /A . l2 

+ [ ^^ + Mi-(i-^] 



—ix(2 — (l — rj)x) 
14 



'i-(i-n)x) 2 



In the narrow gap limit of this equation, for which (1 — rj)x <C 1, we find 

1 /, l-n\ / m 2 (l-r]) 2 \ , , 

-oo xx - U +— gr-^J ( V g " + /i-2u;j£ = 0. (32) 



A. Linear stability limit for small injection angles 

As we saw, instability requires both rotation (the first term of (|32|) ). and accretion (the 
term proportional to e 2 ). Hence we require that the three terms of equation (1321) are of the 
same order of magnitude within the active layer, the layer close to the injection boundary 
(x = 0) where the main development of the unstable mode occurs for small injection angles. 
We bring this out by introducing the scalings 

N 

S = — , = ea and x = ey, (33) 

with N, a and y all 0(1) inside the boundary layer. Equation (1321) then becomes 

j^tiyy -u) y + {a- 2iy)u = 0(e) (34) 

We are interested only in the leading order behavior in e, so we discard the right hand 
side, which is proportional to e. In the narrow-gap limit, the vorticity is linked to the 
streamfunction ip of the perturbation by 

q ~ - Ar + m 2^ = 1 (-tf + m 2 (l - ri) 2 i\ . (35) 



(1 — 77) 2 V e 2 

Hence, at the leading order in e, Co may be replaced in ( )34|) by ip yy . Thus the streamfunction 
close to the boundary satisfies 

jj'&vwv ~ ^wy + ( a ~ 2i v) = °- ( 36 ) 
The solution for ip yy that vanishes for large y can be expressed in terms of the Airy function: 

*m = + Sw + w : ' Nll: 'y) ■ < 37 > 

If we let ipi n denote the corresponding solution that vanishes for large y, i.e., the inner 
solution, we can write the general solution as 

iP(y)=Ay + B + ip in (y). (38) 
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FIG. 5 Inner solution at the onset of instability N = N c = 4.58 and a = 5.62i (solid line: real 
part, dashed line: imaginary part). 

The boundary condition ip(-) = gives Ay + B = A(y — -). The boundary condition 
■0(0) = leads to A = e^i n (0) and, finally, the boundary condition 4> y (0) = gives eipi n (0) + 
(^dy L ^v=° = ®- ^ n ^ ne urm t e ~~ this relation becomes (^*r)j/=o = 0. The parameter a is 
then a solution of the implicit equation 

jf + jg^ + («) * - o- P9) 

The onset of instability is found for N = N c so that a is purely imaginary. We get the 
critical values N c ~ 4.58 and a = 5.62z. The asymptotic behavior for small values of the 
injection angle then corresponds to 

Re c ~ 4.58m 1/2 (l - r?) tan 6~ 3/2 , 
Im(A) = m + 5.62^/mtan(0). (40) 

The imaginary part of the eigenvalue corresponds to the advection of the pattern by the mean 
flow close to the outer boundary: for small injection angles, the eigenmode is a traveling 
wave of dimensionful angular velocity Q. The inner solution at the onset of instability is 
represented in Figure [51 

B. Critical injection angle for the smallest azimuthal wavenumbers 

Close to the critical injection angle for which the critical Reynolds number goes to infinity, 
the instability develops in the whole width of the gap. The parameter e does not necessarily 
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0.2 0.4 0.6 0.8 1 

x 

FIG. 6 Marginally stable mode at the onset of instability e 2 = 0.166 and fj, = 2.04i, for an infinite 
Reynolds number (solid line: real part, dashed line: imaginary part). 



have to be small, but we simplify our calculations by focusing on the lowest of the azimuthal 
wavenumbers — specifically on m 2 (l — r]) 2 <C 1. At the leading order in m 2 (l — 77) 2 , equation 
fl35|) reduces to 

Inserting this into (1321) and taking the limit S — > 00 gives 

e 2 ipxxx + (2ix - /i)i/; xx = (42) 

with solution 

ip xx = exp (-^{fJ.x - ix 2 )^j . (43) 

Imposing the boundary conditions -^((J) = 0, d x (?p) x=0 = 0, ^(1) = 0, we deduce the 
dispersion relation 

/ / exp I -^(/xx — ix 2 ) J dxdt = 0. (44) 

Once again, the critical injection angle is obtained when the solution /1 of this implicit 
equation has vanishing real part. This yields the values e 2 = 0.166 and fi = 2.04i, and hence 
the asymptotic behavior is given by 

tan6 c = 0.166m(l - r]) 2 

Im(A) = m + 2.04m(l -77) (45) 
The unstable mode is shown in Figure [6l 
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FIG. 7 Top: Curves of marginal stability for different values of rj and different azimuthal wavenum- 
bers m, in the e 2 — S plane. The small injection-angle asymptote is valid for every curve in the 
small gap limit. There is always a vertical asymptote but its placement is accurate for small 
wavenumbers only (m 2 (l — r/) 2 << 1). Bottom: (Rescaled) deviation of the angular velocity of 
propagation to that of the outer cylinder. The asymptotic behaviors have been drawn for compar- 
ison. The symbols are o: rj = 0.9, m = 3; □: rj = 0.9, m = 10; o: r/ = 0.9, m = 30; •: rj = 0.99, 
m = 30; +: rj = 0.99, m = 100; •: 77 = 0.99, m = 300. 

C. Collapsing the different marginal stability curves 

The marginal stability curves are plotted on FigureHin the e 2 — S plane for different values 
of rj (in the small gap limit) and different values of the azimuthal wavenumber m. The curves 
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corresponding to the lowest azimuthal wavenumbers (in the sense that m 2 (l — if) 2 <C 1), 
collapse onto a single curve. For small injection angles, this curve admits an asymptote 
given by equation f l40p . Its critical injection angle is given by the vertical asymptote derived 
in (|45[) . For high azimuthal wavenumbers, the asymptote to the limit of small injection 
angles is still acurate, but the vertical one is not. This comes from the fact that one cannot 
simplify the relation uj ~ — (jz^p +m 2 ip ) and thus cannot ignore the dimensionless parameter 
m 2 (l — rf) 2 in the limit of high azimuthal wavenumbers. The critical injection angle goes to 
zero in this limit, and its behavior may be determined with more careful asymptotics. 

The angular velocity of propagation of the unstable modes are collapsed onto a single 
master curve in Figure [7] by plotting the quantity = {^^- — lj /(l — 77) as a function 
of e 2 . This quantity can be thought of as (a rescaled version of) the deviation of the angular 
velocity of propagation to the velocity of the outer boundary. It goes to zero as e tends 
to zero following the scaling law ( |40|) . whereas it tends to the finite value 2.04 as e — > e c 
according to equation ( |45l) . 

V. ENERGY STABILITY OF THE BASIC FLOW 

Energy stability is defined below; it is a strong form of nonlinear stability theory. It 
provides sufficient conditions for a flow to be stable to perturbations of arbitrary amplitude. 
However a flow that is not energy stable may still be linearly stable. In this section we 
formulate the problem of determining the energy stability limit of the base flow introduced 
above. 

A. Evolution equation for the energy of the deviation from the laminar solution 

The starting point for the analysis is the decomposition of the velocity field u into the 
steady laminar flow and a time dependent deviation. In this analysis the deviation from 
the laminar solution need not be small. We thus introduce the decomposition u(r, t) = 
V; am (r) + v(r, t) into equation ([7]) to obtain 

v t + (v. V)v + (V lam . V)v + (v. V) V lam = -Vp + ^Av . (46) 

The time- dependent velocity field v(r, t) is divergence free, satisfies the homogeneous bound- 
ary conditions v = at the two cylinders, and is assumed to be periodic in z with period 
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L z . 

To study the kinetic energy of this variable field we take the dot product of this equation 
with v and integrate over one cell, i.e. over the domain r = [i?i,i?2] X [0, 2tt] x [0,L 2 ]. Let 
us denote the volume element in this domain by dr and introduce the notation 

||f|| 2 = ^|f| 2 dr 

for what is known as the L2 norm of /. Then, after a few integrations by parts that make 
use of the homogeneous boundary conditions on v, we have 

dr = -H{v}. (47) 

If the quadratic form % is strictly positive, i.e., if there is a number [i > such that H{v} > 
fj,\ |v| 1 2 for any divergence- free field satisfying the homogeneous boundary conditions, then the 
kinetic energy of the perturbation v decreases monotonically in time at least exponentially. 
A flow for which H is a positive quadratic form is called energy stable — it possesses a strong 
form of absolute asymptotic nonlinear stability. 
Now we define the function [i(Re, tan 0) by 

fi(Re, tan 6) = inf , (48) 

I l v l I 

the infimum (or greatest lower bound) being taken over divergence-free vector fields satisfying 
the homogeneous boundary conditions. Energy stability is achieved in the region where 
fi(Re, tanO) > in the Re — tan0 plane. The level set /x(i?e,tan9) = is the boundary 
of the parameter region where the flow is marginally energy stable. 



d 
dt 



V 



Re 



Vv| 2 + v.(VV iam ).v 



B. Euler-Lagrange equations 

Only the terms of equation ( I46p which are linear in v contribute to % so that it is clear 
that we can write %{v} = / v • L\dr where L is a symmetric linear operator diagonalizable 
in an orthonormal basis with real eigenvalues. We perform the variational procedure subject 
to the constraint of incompressibility of the flow by adding the integral f pV.vdr to H. A 
second Lagrange multiplier A is used to impose the constraint | |v| | 2 = constant when seeking 
the infimum in equation (j4*8l . We may then seek a stationary condition on the combined 
functional. The Euler-Lagrange equations are then 

Av = -^Av + i[(VV, am ).v + V.(W,am)] + Vp (49) 
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together with the incompressibility constraint V ■ v = 0. As usual, the Lagrange multiplier, 
p, in the Euler-Lagrange equations plays the role of pressure. The Euler-Lagrange equations 
for the variational problem are equivalent to the eigenvalue problem for the linear operator 
L and the infimum in equation (|48p is reached when v is an eigenvector of L for its lowest 
eigenvalue. This lowest eigenvalue is the lowest acceptable value of A in equation (jl9l . We 
stress the fact that although these equations are linear they give sufficient conditions for the 
laminar flow to be stable against perturbations of arbitrarily large amplitude. We solve this 
problem numerically in IVHI but first we consider some analytical aspects. 

In the classical Taylor- Couette flow without suction one can show that for some parameter 
values, the Taylor-vortices — that is a ^-dependent and ^-independent mode — are linearly 
stable but may exhibit transient growth (see for instance Joseph, 1976). In the present 
energy stability analysis of the Taylor-Couette flow with suction, we thus need to retain 
both 9 and z dependent deviation fields v to compute the energy stability limit correctly. 

Since the problem is linear and invariant to translations in 9 and z, it is useful to decom- 
pose v into Fourier modes 

v = (u(r), v{r), w{r))e ime e ikz . (50) 

The modes decouple so we may study the stability of each one separately. Insertion of this 
into the Euler-Lagrange equations yields 



= ( -A + A(r) + t -^) u + Z(r)v+p r - ^—^--(ru r ) r (51) 
\ r z j Re r 

= Z*(r)u + ( -A + A(r) - *^®) v + %- ^ V), (52) 

= ( -A + A(r) - ^=^) w + ikp- {J ^ 11 -(™r)r (53) 



where the functions A and Z are 

l-*7 A a + m2 + 1 
Re V r 



A (r) = tZJL ( k 2 + —3— ) (54) 



and Z* is the complex conjugate of Z. The final equation of the system is the incompress- 
ibility constraint, 

1 v 

-{ru) r + im — h ikw = 0. (56) 
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To determine the energy stability of the flow, we need to solve this system of equations 
together with the homogeneous boundary conditions on v to find the lowest eigenvalue, A. 
The flow is energy stable if A > as can be seen by inserting the modal expression into (H8T) . 



VI. BOUNDS ON THE ENERGY STABILITY LIMIT 

In this section we find some bounds on the location of the curve /x(i?e,tanO) = in the 
Re — tan O plane. To do this, we need not solve the variational problem explicitly, so we 
instead take two other approaches: 

• First we specify a particular choice of initial perturbation and we study its energy 
stability. If the energy of any disturbance can grow, then the flow is not energy stable. 
By identifying such test perturbations we can bound the location of the marginal curve 
of energy stability from one side. 

• Then we derive an analytical lower bound on the quadratic form "H. As long as 
this lower bound remains positive, H is positive and the flow is energy stable. This 
(rigorous) sufficient condition for the flow to be absolutely stable bounds the location 
of the marginal curve of energy stability from the other side. 

A. Stability of the mode m = k = 

Consider variations that are axisymmetric and translationally invariant in the z-direction 
corresponding to the mode m = k = 0. In this case, the mass conservation equation together 
with the homogeneous boundary conditions on u yield u(r) = 0. To find the marginal energy 
stability limit we impose A = in the Euler-Lagrange equations. After multiplication by 
Re /(l — 77) the azimuthal component of the equations simplifies to 



the case a > 1 corresponding to the most unstable situation. We thus obtain the general 
solution. 




(57) 



This equation admits power law solutions 




We are interested in 




(58) 
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The boundary condition v (1) = implies A = and v(rj) = then imposes either B = or, 



more interestingly, ya — 1 ln(r/) = g7r, where g is an integer. As a increases, the first mode 
to leave the domain of energy stability has q = 1 for which the corresponding critical value 

2 

of a is a c = 1 + n^p- ■ O n returning to the Re — tanO plane, we find an upper bound on 
the critical value of the Reynolds number at a given value of tanG, namely 

7T 2 \ 1 — T] 



ReJtan 6) < iWtan 6) = 1 + - . XN0 ^. (59) 

v 1 ~ v ; V ( ln (^)) / tanG v 7 

The curve i?ei(tan(B)) has been drawn on Figure M for several values of the geometrical 
factor r\. This upper bound shows that whatever the injection angle is, energy stability will 
be lost when the Reynolds number becomes high enough. 

Figure |8] illustrates how the m = 0, k = mode may grow — at least transiently — 
when suction spoils the conventional interchange argument for large values of a. A small 
perturbation vq = v(t = 0) in the azimuthal velocity field at an initial radius Rq = R(t = 0) 
will be advected towards the center by the suction. When a is very large, the flow is nearly 
inviscid and the fluid nearly conserves its angular momentum during the inward motion. 
Hence the velocity perturbation becomes 

v{R{t),t) = v -^r, with R(t) < Rq, (60) 

which means that the kinetic energy of this perturbation increases and the mode is not energy 
stable. However, when the perturbation is advected all the way to the boundary layer near 
the inner cylinder, viscous effects dissipate the angular momentum of the perturbation and 
the perturbation is removed by the suction. 




FIG. 8 Expected behavior of an energy unstable mode for m = k = 0: a perturbation at t = (a) 
is advected by the suction and speeds up (b). When it reaches the boundary layer the fluid loses 
its angular momentum and the perturbation is swept away (c). 
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B. Lower bound on the energy stability limit 



Write the quadratic form "H{v} as 



1 — r\ 
Re 



with the indefinite term 



V(r) 



Vv| 2 dr + J^{v} 



tan G 



(61) 



uv — (v — u )dr. 



(62) 



To find a lower bound on TL we observe that v cannot be large within the gap while satisfying 
the homogeneous boundary conditions without also having large gradients. Thus we seek a 
lower bound on J 7 (which may be negative) in terms of the non-negative norm ||Vv|| 2 . 
First of all, for a lower estimate we can drop the positive term ^®. u 2 j n Then using 



the inequality \uv | < |[-« 2 + cv 2 }, valid for any c > 0, we have 



^{v}> 



u 



\r{V/r) 



,, i tan 9 \r(V/r) r 
+ v 2 — — + W 



2. / V r 2 



dr. 



(63) 



Now choose c = c(r) so that the coefficients of u and v are equal: 

X 



c(r) 



where x( r ) = r ^\<y/f) | • This choice implies 



4 

X 2 



(64) 



^{v}> 



tan 



1 + — + 1 ){u 2 + v 2 )dr. 



2r 2 \ y 

We now use the fundamental theorem of calculus and the Schwartz inequality to find 



(65) 



\v r 



I —=V?v r (r)dr 
JRx Vr 



< Win I — 



|tv(r)| 2 fdf. 



The same estimate can be applied to u and yields 



u 

Now we use this estimate in fl65| to find 



+ f 2 <ln^-^-^ J (|w r (r)| 2 + \v r (f) | 2 ) rdf. 



(66) 



^{v} > -||Vv| 



* In ' ' 



Hi 



2r V i? 



1 df. 



(67) 
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Thus "H{v} = ^||Vv|| 2 + J^{v} > if 



V 



ne ^"©GRR^ (68) 

If we compute x( r ) f° r the laminar azimuthal velocity profile found in section III. CI the 
lower bound on the energy stability limit becomes the implicit relation 



a f lln [v) ( L I 



— + 1 df = 1 (69) 



y (r) = 2tan9 , -g- ■ . (70) 



2 J,, r V V " ' X 2 

with 

| ar 2+a _|_ 277 2 + Q | 

The line _Re2(tan B) corresponding to this bound is shown in Figure [9] for several values of 
77. Below this curve, the flow is absolutely stable, so the actual marginal energy stability 
boundary lies above it. 



VII. NUMERICAL COMPUTATION OF THE ENERGY STABILITY LIMIT 

In a study of plane Couette flow with suction, Doering et al. (2000) found that steady 
laminar flow was absolutely stable if the injection angle was above the critical value C ~ 3°. 
At this value of the injection angle the energy stability boundary in the Re — tan G plane 
goes to infinity vertically. In the cylindrical problem, however, the upper bound found in 
IVI.AI clearly rules out such a behavior. It is interesting to see how the energy stability 
boundary evolves from the plane Couette limit 77 — > 1 to a cylindrical geometry with 77 < 1. 
To answer this question, we solved the eigenvalue problem in ( )49l) numerically to determine 
the marginal energy stability curve precisely. 



A. Simplification of the system of equations 



Although the system of equations ( |5T]) - (!56|) shares some similarity with the system of 
equations (!T7l)-(!2~3l studied in section [TTTl the presence of ^-dependence makes its simplifica- 
tion less straightforward. Taking the divergence of the vectorial form of the Euler-Lagrange 
equation and using the incompressibility constraint we obtain 



a 1 / 2/, r , \ \ i m , ,rri \ \ tan© 

Ap = -- (r 2 (V/r) r v) r - — (r (V/r) r ) u + — 



u 

-u r H h 



imv 



(71) 
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Next we apply the Laplacian, A, to (l5Tj) and use the identity 

a / \ /ax Pr 2m 2 
A(p r ) = (Ap) r + F - - —p 



(72) 



to eliminate A(p r ) from the result. The remaining terms involving p are (Ap) r , p r and p, 
which can be expressed in terms of u, v and their derivatives using (ITT]) . floT]) and fl52l) . This 
leads to the first of a system of two differential equations for u and t> , which is 



A 



-2im 1 ,, 2 m 2 + l. 1 / 1 - ?A f 2(1- rj) 

— t. — v + u rr + -u r — (k H - — )u = I — I u rrrr + 



+ 2A + 



/2A 



\ r 



+ 2A r 



Re 
im(V/r) r 



Re r 



A 



+ A„. H - — A[ k + 



2 | m 2 + 1 \_ zm(V7r) rr 2im 



+ 



4im(l — 77) 



i?e r 3 



+ tan6 



f r + 
,2 



/ 6im(l — 77) 



Re r 4 



- Z(k 2 + 



m 



/ \ Re r 2 
2imA\ 
—) V 



k m" \ im im 
~ I T + —T I u + — v r r« 



(73) 



To derive the second equation of the system, we express w in terms of u and v using mass 
conservation ( 1561) . On inserting u> into ( |53l) . we find an expression for p that can be inserted 
into fl52l). This leads to 



A 
+ 
+ 



im 

u r 

r 

im 
r 

1 — 7] 

Re 



'-u+(k 2 + 



im(l — 77) 
Re r 



U r , 



A)ur+\k 2 Z*- t ZA+ 2tm{1 - r]) 



Re r 4 



u 



1 1 m 
k +7T 



'Vrr ~\~ 



(_ k 2 + _____ 

Re r V r 2 



+ A(/r + — 



2(1 -i])m 2 A; 2 tan 6 



+ 



/ 2im(l — if) 
\ Re r 2 



(74) 



This system is of fourth order in u and second order in v and there are, appropriately, four 
boundary conditions on u and two on v. 



B. Prom plane Couette to cylindrical Couette 



We solved the eigenvalue problem numerically using a finite-difference method. The 
resulting marginal energy stability curves are shown for several values of i] in Figure __ 
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FIG. 9 Marginal energy stability boundary for different values of 77 (top: 77 = 0.99, middle: 77 = 0.9, 
bottom: 77 = 0.02). The solid line is the numerically-computed energy stability boundary, the 
dashed line is the upper bound Re\ on the location of that curve, and the dotted line is the lower 
bound Re2- 
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When r] is close to one and the entry angle is very small, the energy stability limit remains 
constant at the energy stability limit of the plane Couette flow, Re ~ 82. When 9 ~ 3° the 
energy stability limit increases greatly but it cannot go to infinity as in the plane geometry. 
When 0.3 < tan G, the most unstable mode is (m = 0, k — 0) and the energy stability limit 
coincides with the upper bound Re 1 . 

When i] goes to zero, the energy stability boundary is very different from that of the plane 
case. This limit corresponds to the situation where the outer radius goes to infinity while 
the inner radius is kept constant, and the boundary curve is a monotonically decreasing 
function of the injection angle. The critical Reynolds number at = is higher in this 
case. 

VIII. CONCLUSION 

It is not unusual to suspect that a given flow may be unstable only to encounter mathe- 
matical arguments to the contrary. A familiar example for fluid dynamicists is provided by 
kinematic dynamo theory. In a working kinematic dynamo, a magnetic field grows exponen- 
tially from a small magnetic perturbation on a suitable flow. To many this prospect once 
seemed to be forbidden by Cowling's antidynamo theorem. The way out of the dilemma was 
to (finally) appreciate that Cowling's theorem posited an axisymmetric flow and to focus 
on flows that did not respect this constraint. In the present problem we confront a similar 
situation: many flows with differential rotation may seem to have enough energy to excite 
the growth of disturbances to the flow yet Rayleigh's criteria appear to forbid a purely fluid 
dynamical instability. In this instance, we have contravened Rayleigh's proscription by sim- 
ply violating one of the conditions of Rayleigh's demonstrations. In both of these examples, 
the value of the relevant (anti-)theorem has been to point the way to making the (seemingly) 
forbidden event happen by vitiating one of the theorem's premises. 

As we mentioned at the outset, there are several interesting problems where differential 
rotation plays an important role. We have not hidden our particular interest in the situation 
of accretion disks where matter rotates differentially around a central gravitating body. 
These disks are often quite luminous and the source of their emissions has been thought to 
be in the dissipation of the turbulence that they had long been assumed to sustain. But 
these objects typically have Keplerian rotation laws where their circular velocities vary as 
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1/y/r, where r is the distance from the central object. Hence there was no generally accepted 
purely fluid dynamical process to produce this turbulence, though the matter was steeped 
in controversy. Moreover, turbulence in the disk, if it occurred, was supposed to transport 
angular momentum out of the disk. This would allow the purely circular motion to begin 
to turn inward so that the gas in the disk could fall onto the central object and be accreted. 
The resulting spiral motion, or swirl, is then a central feature of the accretion disks yet its 
dynamical role has been largely glossed over in the discussions of disk instabilities. 

For us, the intriguing feature of the instability of swirling (a.k.a. spiraling) flows in disks 
is that their instabilities may promote accretion that may in turn enable the inward flow 
that was the source of instability. This is positive feedback par excellence and we may 
reasonably anticipate that the instability treated here is likely to cause subcritical bifurcation 
in disks, though we have not yet shown that. Nevertheless, this physical argument leads us to 
suppose that rather than thinking of a competition among instability mechanisms, we may 
imagine a cooperation amongst them. Thus in the linear study of the magneto-rotational 
instability Kersale et at. (2004) find that the resulting fluid stresses drive accretion for 
suitable boundary conditions. This self-induced accretion means that a swirl is produced 
that is likely to contribute to the instability we consider here. This interaction of instabilities 
in the fully nonlinear regime bids fair to give rise to many of the processes that have been 
sought to explain the observed wealth of behavior in swirling flows. 

However, the work reported here does not treat disks explicitly. We have chosen instead 
to consider a fluid flow field that is well known to be stable and to show that by converting 
it from purely circular to swirling motion by the agency of accretion we could render it 
linearly unstable. This required only a simple change in the boundary conditions to vitiate 
the applicability of Rayleigh's inflection point theorem for circular motion. It remains to 
learn whether this destabilization will apply to Keplerian flows, though the mechanism seems 
generic and we do think the instability is inherent to physically realistic, spiraling flows. As 
to real disks, there we are in the domain of shallow gas theory and that will lead to a more 
difficult investigation, but in the beginning it will be a linear one. In the meantime, it 
could well be that other instabilities will be uncovered as other variations on disk flows are 
examined. These will also join the magnetorotational instability that now leads the pack 
and whose consequences are currently being vigorously explored by many investigators. 
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Appendix A: Bounds on the energy dissipation 

The first physically relevant rigorous bound for turbulent flow was given by Howard (1963) 
using an approach to thermal convection proposed by Malkus (1960). Subsequently, related 
techniques have been developed such as the background method (Doering and Constantin, 
1994, 1996). Although the background method can be applied to many problems in many 
geometries, there has been much emphasis on shear layers. For example Constantin (1994) 
used the method to compute a bound on the energy dissipation for the Taylor- Couette 
flow without suction and Doering, Spiegel & Worthing (2000) found a bound on the energy 
dissipation in a plane shear layer with suction. However, at high Reynolds numbers and 
strong suction the Taylor-Couette problem appears to be significantly more challenging. We 
show here that the usual formulation of the background method is inadequate in this case. 

1. The background method 

The background method begins with a decomposition of the velocity field into a steady 
velocity profile and a time-dependant fluctuation field u(r, t) = V(r)+v(r, t). The difference 
from energy stability analysis is that the background profile remains to be chosen; it is not 
necessarily the basic steady solution. Moreover it needs not be a mean field as in the 
Malkus-Howard approach. We insert this decomposition into equation (JTj), take the dot 
product with v and integrate over the volume r. We perform some integrations by parts 
and use the equality 

|Vu| 2 = | VV| 2 + | Vv| 2 + 2(Vv) • (VV) (Al) 

to deduce 
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with 



1 — 7] 
2Re 



Vv| 2 + v.(VV).v + v.(V.V).V 



dr. 



(A3) 



The last term in I is linear in v which is not desirable in this procedure so, to alleviate this 
problem, we introduce another decomposition, 



v = W(r) + w(r,t) 



(A4) 



where W and w are both divergence-free and satsify the homogeneous boundary conditions. 
(W is steady whereas w may be time dependent.) This decomposition is inserted in / to 
yield 



/ = X{w} + ^^||VW|| 2 + / [W- (VV) • W + W- (V- V) • V]dr + £{w} (A5) 



2Re 



where 



X{w} 
£{w} 



l — i] 
2Re 



Vw| 2 + w- (VV) • w 



dr 



w 



1 — 7] 

Re 



AW + (VV) • W + (W • V)V + (V • V)V]dr. (A6) 



If the quadratic form X is positive, a field W exists that may be employed to remove the 
linear term in w (X > is an invertability condition for the linear operator acting on W 
inside the square brackets in £). Then £{w} = for any vector field w that is divergence- 
free and satisfies the homogeneous boundary conditions (this does not mean that the bracket 
inside C has to be zero, just that it is a gradient). When such a vector field W is found, 
then £{W} = as well, and 



J{W} = ~ J W • (V • V)Vdr. 

The quantity / is then 

/ = X{w} + i J W • ((V • V)V)d 
which when inserted in equation ( 1A2I) yields 



(A7) 



(A8) 



d t (||v|| 2 ) + i^HVul 12 - 1 



Re 



Re 



VV| | 2 - / W • ((V ■ V)V)dr - 2X{w}. (A9) 



31 



X{w} is the same quadratic form as for the energy stability analysis of V as if it were a 
steady solution with Re replaced by 2Re. If the background velocity field V may be chosen 
so that X is a positive quadratic form, then the corresponding W exists and X{w} can be 
dropped in (IA9I) so that 



where the overline represents a long time average. 

In most applications, background profiles are chosen to depend on the same coordinates 
as the simplest steady laminar solutions that respect the symmetries of the problem. By 
choosing the background as a function of the depth only, we may concentrate the background 
shear into thin layers near the no-slip boundaries. Background flows can then be found that 
ensure that X is a positive quadratic form, at the cost of increasing the bound on the right 
hand side of ( lAlOj) . Choosing an azimuthal background velocity profile that depends only on 
r as for the laminar solution is a successful strategy for the classical Taylor Couette problem 
(Constantin, 1994). 

In the Taylor-Couette flow with suction, there is a limit on the energy stability domain of 
any background flow that cannot be extended: the energy stability boundary that emerges 
from the study of the (m = 0, k — 0) mode, does not involve the azimuthal velocity profile. 
This means that for any injection angle B, and for any Reynolds number greater than 
2Rei(Q,rj), the quadratic form X will not be a positive quadratic form no matter what 
function of r is chosen to be the background azimuthal velocity profile. Therefore we cannot 
produce a bound on the energy dissipation with a "simple" background azimuthal velocity 
depending only on r. (This sort of mathematical obstruction was also recently discovered 
in an abstract setting by Kozono and Yanagisawa (2009).) There is no doubt a message in 
this limitation on the method but we shall not speculate on this here. 

In studies on convection or shear layers, optimized background profiles resemble mean 
flow profiles at high Rayleigh and Reynolds numbers (Doering and Constantin, 1994 and 
1996). This reinforces the common belief that the symmetries of the problem which are lost 
through symmetry-breaking bifurcations are recovered in a statistical sense for developed 
turbulence at even higher Reynolds numbers. The fact that we cannot choose a background 
velocity profile for the Taylor-Couette flow with suction that respects the symmetries of the 
problem raises the question of whether or not these symmetries may be recovered in any 




(AlO) 
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sense at high Reynolds number. 
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